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Abstract 

We investigate the dynamics of matter and optical waves at the early stage of superradiant 
Rayleigh scattering from Bose-Einstein Condensates. Our analysis is within a spatially dependent 
quantum model which is capable of providing analytic solutions for the operators of interest. The 
predictions of the present model are compared to the predictions of a closely related mean field 
model, and we provide a procedure that allows one to calculate quantum expectation values by 
averaging over semiclassical solutions. The coherence properties of the outgoing scattered light are 
also analyzed, and it is shown that the corresponding correlation functions may provide detailed 
information about the internal dynamics of the system. 

PACS numbers: 03.75.Nt,67.85.-d,37.10.Vz,42.50.Ct 
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I. INTRODUCTION 



Superradiance in general terms refers to enhanced emission from an ensemble of radiators. 
It was first predicted by Dicke in 1954 [l| and since then experimentally confirmed to occur 
in many systems, such as gases of excited atoms, molecules or quantum dots 0, ^ Re- 
cently, superradiant scattering off an elongated atomic Bose-Einstein condensate (BEC) has 
received much theoretical (4 - T ] and experimental 8-1(3] attention. There are many similar- 
ities, but also important differences between the "conventional" superradiance, for example 
off excited gases, and the superradiance off atomic condensates. Atoms in a BEC have a 
narrow momentum distribution, and thus the recoil they experience during the absorption 
and emission of photons has a profound impact on their momentum distribution, leading to 
distinct atomic scattering patterns. 

In the case of superradiant Rayleigh scattering off BECs, different regimes of parameters 
have been identified, which are characterized by distinct atomic patterns. Mean-field models 
were found to successfully predict and explain such patterns, as well as the transition between 
different regimes, provided the models include spatial effects 4j, |6[. The main drawback of 
such models, however, is that one has to "seed" the equations of motion, in order to start 
their evolution in time. The seeding introduces some ambiguity in the solutions, which is 
expected to become less important for large times due to the fast growth of the population 
in the various optical and mater-wave modes. In contrast to the mean-field models, the 
quantum models that have been used in this context, are capable of describing accurately 
the startup of the process, but do not take into account spatial propagation effects 5, ij]. 



In a recent work 



11| , we investigated the coherence properties of matter waves produced 



in superradiant scattering off BECs, and analyzed the type of spatial correlations involved. 
This has been possible in the framework of a spatially-dependent quantum model, which can 
describe quantum fluctuations while capturing spatial effects, essential for a full understand- 
ing of the process. The purpose of the present paper is to provide a full derivation of the 
model used in H[, and to obtain further detailed insights into the dynamics of the system. 
We explicitly compare our results to those obtained within a related mean-field model, and 
show that for a large collection of condensed atoms, the effect of quantum fluctuations on 
various observable quantities can be obtained by averaging over ensembles of semiclassical 
solutions ("trajectories"). Finally, we discuss the temporal coherence of the scattered light, 
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which is shown to carry information about the internal dynamics of the system. 



II. THE MODEL 

The system under consideration pertains to a BEC elongated along the z-axis, consisting 
of N atoms. A linearly polarized laser with frequency ui = kic, far detuned from the closest 
atomic transition by a value of 5, is illuminating the cloud along the x-axis. Assuming two- 
level atoms and adiabatically e . 
equations for the system read 



iminating the excited atomic state, the Maxwell- Schrodinger 



K5 



(lb) 



with the atomic polarization 



pH( Xl t) = -d^( x ,t) d " E( ^ (x ^^ (x,t), (lc) 

satisfying p(+) = p(~' . The operators E^^x, t) are the positive and negative frequency 
parts of the electromagnetic field, while ip{?t,t) is the operator describing the ground state 
of the atoms, which have mass M and dipole moment d. In writing Eqs. (JT|) we have 
neglected the external trapping potential as well as atomic interactions, which both do 
not play a significant role for the timescales of interest. Due to the coherent nature of 
the condensate, successive Rayleigh scattering events are strongly correlated and lead to 
collective superradiant behavior. As a result of the cigar-shape of the condensate, the gain 
is largest when the scattered photons leave the condensate along its long axis, traveling in 
the so called endfire modes with wave vectors k « ±kie z . A condensate atom can scatter a 
laser photon into the endfire modes, experiencing a recoil h(\ ~ fr(kie x — k). On the other 
hand, it can also scatter a photon from the endfire modes into the laser mode, in which 
case its momentum changes by hq ~ h(—kie x + k). These processes lead to the formation 
of two pairs of atomic side-modes, consisting of counterpropagating atoms with a narrow 
momentum spread (compared to k{). Of course, atoms within these side- modes can also 
scatter photons, thereby acquiring higher momenta, but since we are interested in the early 
stage of the process, we consider only first order sidemodes to be populated. 



Neglecting any coupling between counterpropagating photonic endfire modes, the system 
becomes symmetric with respect to the rr-axis. We can thus focus on the endfire modes with 
k ^ +ke z and on the two atomic side- modes (with central momenta hq = ±h(kie x — ke z )) 
that are coupled to them. Due to the strong confinement along the x- and y-axis, we can 
assume the transverse profiles of the matter field ip± (x, y) to be well described by a classical 
function, independent of the ^-coordinate 6| . Assuming a Fresnel number close to unity for 
the electromagnetic fields, we can apply the same approximation for the transverse part of 
the radiation field u±(x,y), effectively reducing the problem to one dimension. 

We expand the field operators as 

1 

$(x,t) = 4)±{x,y) 4-M)e y(fc,a, ~ fc * Hw '* (2a) 

E (+) (x,i) = ^e y e^-^ +u ± (x,y)E ( +\z,t)e y , (2b) 

where u±i = H(kf + k 2 )/2M and uq = 0. The matter- wave operator is split up in three 
parts, describing the two side-modes (j = ±1) and the BEC at rest (j = 0), 

^(z, t) = e 1 ^ -7=c. jk+p (t), (3a) 

P 6Ao VL 

where A is the interval (— k/2, k/2) in /c-space, L is the length of the BEC and c p annihilates 
an atom with momentum Hp. Since the BEC at rest remains practically undepleted it can 
be treated as a time independent classical function and hence we can set ipo(z,t) = ipo(z). 
Similarly, we expand the endfire mode operator as 



£«(z,t) = ie^-*) £ J^^e^a k+P (t) 

peAo v 

~ l^e^—t)e + ( Zyt ), (3b) 
V 2^0 

where a p is the photon annihilation operator. The frequencies w&+p are approximated by 
u = k/c, the frequency of the scattered photons, and can therefore be taken out of the sum. 
This approximation is justified by the fact that dominant contributions to the sum come 
from momenta of order 1/L, which is several orders of magnitude smaller than k. 
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The one- dimensional field operators satisfy the commutation relations 



e + (z,t 1 ),eUz,t 2 ) 



5 ij 5 A (z 1 - z 2 ), 
- z 2 ), 
#a(*i - 1 2 ), 



'IJ 

5a(^i 

1 

c 



(4a) 
(4b) 
(4c) 



where 8&(z) is a distribution with width of order 1/k and Sij denotes the Kronecker delta 
Inserting expansions (j2D in the equations of motion (CQ), we can make some further 
simplifications. Considering that ki is about a factor of 10 smaller than the extent of the 
BEC cloud along the strongly confined axes, the two transverse functions ip±(x,y)e ±lklX are 
mutually orthogonal to a very good degree of approximation. Hence we can project on 
either of the two side-mode operators by multiplying with the complex conjugate of the 
corresponding transverse function and integrate over the variables x and y. 

Since we included the phase factors arising from the free time evolution in the definition 
of the operators, we can apply the slowly- varying-envelope approximation (SVEA), which 
yields the equations 



dr 

0^-1 (£,t) 



dr 



+ X 



dr 
<9e + (£,r) 



at 



i/se+(6> T )^o(0> 
-i«e + (£,T)V>o(0-2i^_i(£,T), 
-i [Wo(0$.i(£,t) 

+«#-i(£,t)Vo*(0" • 



(5a) 
(5b) 

(5c) 



Here we have defined t) = e ?/>_i(£, r) and rescaled length and time to dimension- 



less units 



£ = fcjZ, r = 2u; r i, 



(6) 



where u r = Hkf/2M. Accordingly, the fields are rescaled as 



e+(£,r) 



e + (z,t) ipj(C,r) 



ipj(z,t), 



and the speed of light becomes x = The effective one- dimensional coupling is given by 
k = gy/kiL/ (2u r ) with 

d ■ e y \ 2 t 
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h 2 5 



fioj 



dxdyu ± (x,y)^ ± (x,y). 



The SVEA pertained to neglecting derivatives of slowly varying functions in order to arrive 
at Eqs. (]5a]) . fl5b|) . For Eq. fl5c|) . we have kept first order derivatives, but neglected terms 
proportional to e + (z, t) in comparison to the laser field. The frequencies in our system satisfy 
u±i <C 0Ji, bj and thus we only kept time derivatives involving e~ luJlt and approximated u ~ ui. 

Backwards recoiling atoms are a particular feature of superradiant Rayleigh scattering 
off condensates. The physical process underlying the backwards modes violates energy con- 
servation by an amount AE ~ 4ftw r . Thus, according to Heisenberg uncertainty principle, 



it can take place only for times shorter than a critical time t c = h/ AE [9], which, in our 
units, is given by 

r c = 0.5. (7) 

For such short pulses, and for sufficiently high power one typically observes an X-shaped 
pattern for the distribution of the atomic side modes with the initial BEC in the center 
and the recoiling atoms moving both in and against the direction of the applied laser pulse 
(Kapitza-Dirac or strong-pulse regime). On the other hand, for weaker pulses with duration 
longer than r c , the distribution of the side modes exhibits a fan pattern, involving mainly 
forward recoiling atoms (Bragg or weak pulse regime). If we neglect the atomic backwards 
sidemode altogether, Eqs. ([5]) become formally equivalent to descriptions of "conventional" 
superradiance from excited atomic gases jlj]. 

Finally, the equations of motion fl5]) can be derived from the effective, self-adjoint Hamil- 
tonian 

A = I ^ + + /#0f ^ 1 + "^W-i + hx -) ' (8) 

where "h.c." stands for the Hermitian conjugate. The system being effectively hamiltonian 
guarantees conservation of the commutation relations (j4j) for all times. Differentiating the 
atomic densities 

n j {Z,T) = (ty{Z,T)$ j {£,T)) (9a) 

and the photon density 

a(£,r) = (ei(£,T)e + (£,r)> (9b) 

with respect to time, using Eqs. ([5]) and adding up the resulting three equations, we find 
the continuity equation 

— [n +1 (£, r) - n_!(e, r) - 3(£, r)] = r). (10) 



Let us now integrate this equation over time from to r and over space from one end of the 
condensate at £ = to the other end at £ = A = kiL. Assuming that side-mode and photon 
populations vanish at r = 0, we find 



Af +1 (r) - A/li(r) - 2jn(r) = X out (r) 



(11) 



where we have defined the total populations for the atoms 




(12a) 



o 



and the photons 




Xout(r) = X / dr'3(A,T'). 



d£3(£,T) 



(12b) 



(12c) 



In words, Eq. (Hip expresses that at any time r, the number of forward-recoiling atoms 
A/+(r), is equal to the sum of backwards recoiling atoms A/l(r), and endfire photons inside 
and outside the BEC volume, denoted by X in ( out )(r). It is therefore consistent with the intu- 
itive picture of the underlying process and it may serve as a convenient check for numerical 
simulations. 

III. SOLUTIONS OF THE EQUATIONS OF MOTION 

We can use the Laplace transform to find exact solutions to the system (jSJ) in terms of 
the operators evaluated at the boundary of their domain - i.e. at £ = and r > or vice 
versa at £ > and r = 0. More details on this procedure are given in the appendix. The 
solutions read 
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IK 



dr'e + (0, r')F 0fi (^ , r - r' - - - / Vo(0^+i(0 O^ofe, r - 
o X Jo 



£ 



+^o*(0^-i(0 0)^0,1(7^', r - %0 + -e+(0 0)F ,o(7«', r - fa 



K 



(13a) 



^i(e,r)=i«^(0 / rfr'e + (0,r / )F li0 (7 ? , ,r-r / -/3 5i o)+^ 1 (e,0) 
Jo 

-2 K 



+ -^o(0 / ^0(0^1(^,0)^,0(7^^-^) 
X Jo L 

+ v> *(O^-i(O o^^, r - & j£ o +-mo 0)^1,0(7^, ^ - 

i>U(t, r) = - i«Vo(0 T ^e + (0, r / )F , 1 ( 7€ ,o, r - r' - &, ) + e^L^, 0) 

Jo 

- -^o(e) / k(0$-i(0 o^ife, r - p u .) 

X Jo 1 



(13b) 



+^(O^(O0)F 0>2 (7^,t - %0 + -e + (e',0)F 0il (7^,r - fa,) 



(13c) 



where we have introduced 



7^ = -[p(0-p(O], 



with p(0 = lodC'IMC )| 2 - The functions F^u, u) are defined as 



^(w,^) = £ * 



e u/p e -u/(p+2\) 



where C~\ v denotes the inverse Laplace transform. One can check easily that Eqs. (I13p 
indeed are solutions to the system (jSJ) by using recursion relations for the functions F^ v which 
are given in the appendix, alongside the explicit expressions of the functions themselves. 
Explicit expressions for the functions F^ u {u,v) appearing in ( 1T3|) are given in the appendix. 
They are combinations and integrals over combinations of Bessel Functions. It shall only 
be noted here, that all of the terms appearing in F^ u {u,v) contain Heaviside step functions 
0(f), except one term in i*o,o(w, v), which is simply the Dirac delta function with argument 



v. 



We note that in Eqs. ( TTBl) . all time arguments are shifted by the value /%£', which is the 
time a photon needs to travel from £' to £. Using the step functions in the solutions to change 
the range of the integrals and assuming free light propagation outside the condensate we can 
reformulate Eqs. fflBl . such that they involve only spatial integrals ranging from £ — t\ to £. 
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This is a consequence of the finite speed of light, allowing atoms at £ only to be influenced by 
atoms within a range £ — r%. Such retardation effects are very small in the system at hand 
and can be neglected for all practical purposes. We can do so formally by letting x — > oo, 
which implies P^^> — > 0, and neglecting all the terms proportional to e + (£', 0) in the spatial 
integrals of Eqs. (fT3"j) . The resulting solutions will still describe the system correctly, since 
the effects of this approximation are expected to be of the order of A/% ~ 10 -7 , and are 
thus too small to be noted in typical BEC experiments. Formally, this approximation will 
lead to a nonzero initial photon population within the BEC, which we can safely neglect 
due to its small value. 

Equations (1131) can be simplified considerably if we neglect backward recoiling atoms. 
Neglecting retardation effects, we obtain 

e+(£,r) = / dr'e + (0,r')Fo(^,o,T-r') 
Jo 

-- /VlMO^-iC?, °) F ^',r) (14a) 
X Jo 

^(£,r) = h#*(£) [ T dr'e + (0,T')F 1 (^ ,T-r') + ^ +1 (^0) 
Jo 

+ -r o (0 [ «o(0^i(e',0)F 2 ( 7 ^,r), (14b) 
X Jo 



where F^(u,v) = C~\ v |e"/ p p~ M }; explicit formulas for F^iu.v) are given in the appendix. 
It is worth emphasizing that these equations are consistent with the equations other authors 
derived to describe conventional superradiance 13]. 

Assuming that the initial population of the atomic side modes is zero, we can find the 
expectation value of any correlation function pertaining to electromagnetic- or matter-wave 
fields using Eqs. (j4|) and ( Ti~3i) and calculating the occurring integrals numerically. For 
the numerical calculations, we assumed the BEC to consist oi N = 10 6 Thomas-Fermi 
distributed 87 Rb atoms, such that ^ = y^Q(z)NQ(Lz — z 2 )/L 3 with L = 130/im. We used 
a spatial grid of 400 points. For the incoming laser we chose a rectangular profile and a 
wavenumber k\ = 8.05 x 10 6 m~ 1 , which results in a dimensionless length of the BEC of 
A ~ 1000. Coupling strengths are conveniently expressed in terms of the superradiant gain 
T = k 2 N/x, whose value separates the two regimes identified by experimental observations 
of superradiance from condensates [h| . Typically, the weak coupling regime is characterized 
by g ~ 10 5 s _1 and T ~ 1, while for g ~ 10 6 s _1 and T ^> 1 the system is in the strong coupling 
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regime. In our calculations, we chose T = 1 and T = 100 for the two regimes. 



IV. QUANTUM VS. MEAN FIELD DESCRIPTION 

Various aspects of the strong and weak-coupling regimes have been described successfully 



within a mean- field (MF) model jg], which is closely related to the present quantum model 
given by Eqs. (jSJ). In fact, one arrives at Eqs. (15) of js] by adapting the approximations of 
the present quantum model and replacing the operators in Eqs. §5§ with their expectation 
values, treating them as classical fields. Consequently, we will refer to the solutions of the 
mean field model as r) = t)) and e+(£, r) = (e+(£, r)). Due to the generality of 

the Laplace transform, these solutions look exactly like Eqs. (|T3|) . but with the operators 
replaced by classical fields. 

Both models take into account spatial effects, which have been shown to play a major role 
in Rayleigh superradiance from condensates Given, however, that our system is initially 
prepared in the vacuum state, both ipj(£,,r) and e+(£, r) will remain zero throughout the 
evolution of the system, because their equations of motion [see Eqs. (|T3|) ]. are not coupled 
to any operator with non-zero expectation value [i.e., 0) = and e+(£, 0) =0]. This is 
a major drawback of the MF model, which can be resolved by seeding either of the modes 
i.e., assigning a non-zero initial value to either ?/^(£,0) or e+(£, 0). The arbitrariness of 
such a seeding introduces some ambiguity regarding the dynamics of the system for short 
times, where all the modes are scarcely populated. The MF model is expected to be valid for 
longer times, where the fast growth of the population in the modes eliminates any ambiguity 
caused by the initial seeding. At such times, the MF model explains reasonably well various 
experimental observations (6j. Our purpose in this section is to investigate how accurately 
one can describe the onset of superradiance from condensates, in the framework of this 
MF model. To this end, we will compare the predictions of the MF model for various 
observables, to the corresponding predictions of the quantum model, which is capable of 
describing accurately initial quantum mechanical fluctuations, and does not require any 
seeding. 
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A. Atomic Densities and Populations 



A rather convenient observable in a superradiant scattering process from a BEC, is the 
atomic density of the sidemode j, denoted by nj(£,r). After turning off the atomic trap, 
atoms in the sidemodes separate from the BEC at rest due to their additional momentum 
and form observable scattering patterns 9|. 

In the quantum model, we have rij(£, r) = (?/>](£, r)V>j(£, r)), which in view of Eqs. f|T3|) 
yields 

n-i&r) =PV(0| 2 /Vk(Ori^,i(7^,r)| 2 , (15a) 
Jo 

n+i^r) =n_ 1 (^r) + T\ V (0\ 2 f dr'|F li0 ( 75 ,o, r')| 2 , (15b) 

Jo 



where = ipo(£)/vN. In the MF model, expectation values are defined as the squared 
modulus of the classical functions, i.e. nj(£,r) = t)| 2 . By means of these quantities, 

we can directly compare the two models. 

As a first step, let us neglect for the time being backward recoiling atoms. In this case one 
can obtain analytic expressions for the atomic densities, which acquire particularly simple 
forms for a fiat BEC [i.e., for ipo(0 = y/N/A]. Equation ( 115bj) reduces to 



n +1 (£, r) = Tr [/ 2 (2 y/T&) - i?(2^^)J , (16) 
with r = T/A, whereas within the MF model one obtains 

n+l({ir) = l*«M!! J j (2 ^ ) . (17) 

Note here that the prefactor of the Bessel functions in the case of the quantum model is a 
linear function of time, as opposed to the time-independent variable 0)| 2 /A in the 

case of the MF model. This is a key difference, whose implications become clearer if we look 
at the asymptotic behavior of the total population of the forward atomic side mode. Using 
Eq. (EE6D in Eq. ( fl2al) . we obtain for Tr > 1 

A/" + i(r) ~ TTT^e 4 ^, (18) 

l07TV 1 T 

while for the MF model one finds 
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where we have assumed a spatially independent seeding i.e., 0) = y/fj. Clearly, 

as a result of spatial propagation effects, both models predict a sub-exponential growth 
of the side-mode population. In the quantum model, however, the population grows like 
e^^/y/r, whereas in the MF model it increases as e 4v/7 /r. The crucial point is that we 
cannot compensate for such a difference by assigning any constant value to the seeding 
77. Furthermore, the fact that the seeding appears as a prefactor in Eq. ( fl6l) . suggests 
that any deviations of the MF atomic density profiles and populations from their quantum 
counterparts have to be attributed to the semiclassical nature of the MF model and not 
to the arbitrariness of the initial seeding, which may only lead to global changes such as a 
rescaling of the plotted curves. Keeping this in mind, we turn to comparing the predictions 
of the two models taking into account both forward and backward recoiling atoms, as well 



as a Thomas- Fermi distributed BEC. For direct comparison to previous work jg], we have 
decided to seed the mean-field model according to 



^ +1 (£,0)=Vo(0/ViV, (20) 

which corresponds to one atom in the forward atomic sidemode distributed proportionally 
to the density of the BEC. 

A snapshot of the atomic density profiles in the two models, after a time Tr = 6, is 
plotted in Fig. [TJ Both models predict that the profiles are peaked close to the right end of 
the condensate; a feature which is responsible for the experimentally observed asymmetry 
of the scattering pattern 6|. In the MF model, however, the profiles are peaked slightly 
closer to the end of the BEC, while the height and the width of the spatial distributions 
are underestimated, especially for the (— ) mode in the weak pulse regime. This discrepancy 
in the predictions of the two models is also reflected in the time evolution of the atomic 
populations in the two side modes, which are depicted in Fig. [21 

In agreement with experimental observations, both models predict that the populations 
of the two side modes are comparable in the strong-pulse regime, whereas in the weak- 
pulse regime, we have far less backwards than forwards recoiling atoms. This behavior is 
mathematically mirrored in the expressions for the atomic densities in the two side modes 
[see Eq. (1151) ]. They differ by one term only, which is proportional to T, whereas their 
common term scales with T 2 . Thus, for short times and strong pulses where T ~ 10 2 , the 
two expressions become comparable, whereas they are different in the weak pulse regime 
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Weak Pulse (r = 1) Strong Pulse (r = 100) 




£ 10 £ 1 

A A 



FIG. 1: (Color online) Atomic density profiles of the two side modes, according to the quantum 
model (black, solid) and the MF model (red, dashed), at Tt = 6. The left column shows the weak 
pulse regime and the right column the strong pulse regime. 
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FIG. 2: (Color online) Evolution of the side-mode populations in the quantum (black, solid) and 
the MF (red, dashed) model. Note that the quantities are plotted on a log-scale. The dot-dashed 
line marks the time when the snapshots in Fig. [I] were taken. 

where T ~ 1. It is also worth pointing out here that, as depicted in Fig. [5J the MF model 
gives approximately the right growth rates as well as the right qualitative behavior. In the 
quantum model, however, the suppression of the (— ) mode is not as prominent as in the 
MF description, and this can be attributed to the ambiguity of the seeding and therefore 
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FIG. 3: (Color online) Comparison of the predictions for the photon density within the BEC in 
the quantum (black, solid) and the MF (red, dashes) model. The snapshots correspond to Tt = 6, 
in the weak-pulse (r = 1) and the strong-pulse (r = 100) regimes. 

the initialization of the process. 
B. Scattered light 

More insight into the differences between the two models is obtained by also studying the 
behavior of the radiation field. In the quantum model, the photon density within the BEC 
volume model is given by Eq. ( I9bl) . which in view of Eq. (I13a,[) yields 



Even though this quantity cannot be measured, we will use it to compare the two models as it 
influences many measurable features of the process. The most easily measurable observable 
linked to the radiation field is the number of photons which have left the condensate up to 
time r. Assuming no distracting factors between BEC and detector as well as instantaneous 
photon propagation, this quantity is given by Eq. (I12cp . Finally, it is straightforward to 
find the analogous quantities in the MF model using 3(£, r) = |e + (£, r)| 2 , and again we can 
directly compare the predictions of the two models. 

Calculations of the photon density within the BEC are plotted in Fig. [3j In the weak 
pulse regime, the mean-field model shows reasonable qualitative agreement with the quantum 




(21) 
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FIG. 4: (Color online) Comparison of predictions for the number of photons which have left the 
BEC as a function of scaled time 1~Y in the quantum (black, solid) and the MF (red, dashed) model. 
Top panel shows the weak-pulse regime (r = 1), while the lower panel is for the strong-pulse regime 
(r = 100). The dot-dashed line marks the time when snapshots in Fig. are taken. 

predictions. In the strong pulse regime, however, the quantum model predicts a much lower 
photon density than the MF model. This discrepancy becomes even more obvious if we 
look at the number of emitted photons, which is plotted in Fig. HI While the MF model 
predicts a fast growth of the number of photons, the growth in the quantum model is almost 
linear, which in view of Eq. (\12c\i implies a constant density of photons within the BEC 
volume. Upon closer investigation, we find that such a period is also present in the weak 
pulse regime, albeit for shorter (scaled) times. More precisely, we find that for T = 1 this 
period lasts only until about r « 0.5, as can be seen in the inset of fig H] (a). 

According to our simulations, the presence of backwards recoiling atoms is suppressing 
superradiance. As depicted in Fig. HJ the number of scattered photons with respect to 
scaled time in the weak pulse regime is much higher than in the strong pulse regime, since in 
the latter endfire photons are destroyed on account of producing backwards recoiling atoms. 
This removal inhibits the fast growth of the endfire mode, which in turn is responsible for 
the lower scattering rate (per scaled time). In particular, the endfire mode remains weakly 
populated for times r < r c , since in this regime, the production of backwards recoiling atoms 
is allowed. This behavior is consistent with the conservation law (llljl . which says that the 
number of photons outside the condensate is given by the number difference between the two 
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FIG. 5: (Color online) Example of photon and atom densities obtained through averaging over 
randomly seeded MF solutions. Figure (a) shows a snapshot of the photon density and figure (b) 
shows the atomic density of the (—1) mode. The quantum solution is shown in solid black, while 
averages over 20 and 2000 trajectories are shown in red dashed and blue dotted lines. Parameters: 
r = 1, and t = 2. 

matter-wave modes. Finally we note that the suppression of the population of the endfire 
mode seems to be underestimated in the MF model, as can be seen from Fig. H] (b) as well 
as the inset of Fig. @] (a). 



C. Averaging over semiclassical trajectories 

For early times, the quantum prediction of the superradiant process is appropriate and 
can be expected to give better results than the MF model. It is easier, however, to perform 
calculations involving depletion of the BEC and population of higher modes in the MF 
model. It is therefore reasonable to ask, whether the quantum model is able to give us 
hints on how to seed the MF model appropriately to obtain quantitatively correct results. 
In the context of conventional superradiance, Haake et al. introduced the idea of averaging 
over many semiclassical "trajectories" to obtain quantum results 13]. The MF equations 



are initially seeded with random variables according to a particular distribution, while to 
obtain a particular quantum expectation value, one has to average over various solutions for 
the corresponding semiclassical quantity. We have investigated the extension of this idea to 



1(3 



superradiant Rayleigh scattering off BECs where, in contrast to conventional superradiance, 
backwards recoiling atoms are also present. 

To see how this works, let us assume we want to calculate the normal-ordered rath order 
correlation function 

(^1i(£i,t)J [V-i(6,t)J }• (22) 

From Eqs. (fT3j) and using fll]) as well as the fact that our initial state is the vacuum for all 
modes, we find the only non-vanishing expectation value involved to be 

0) • • • ^ + i(£ (n) , O)^^, 0) . . . ^(e (2n) , 0)). (23) 

Here, the variables are integrated from to £i for j = 1, . . . ,n and from to £2 for 
j — n + 1, . . . , 2n. Using the commutation relations (jl]), correlation ( 1231) reads 

n 

En^ a) - e (n+7r °' )) ), (24) 

n j = l 

where 5 denotes the Dirac delta and the sum runs over all permutations it of order n. Let 
us now seed the semiclassical model with 

?M£,0)=Q, (25a) 
V_i(£,0)=0, (25b) 
e+(£,0) = 0, (25c) 

where is a random, normally distributed complex variable, with zero mean and variance 
l/yA£, with the spatial step of a numerical implementation. The average of the prod- 
uct ?/>+i(£ (1) , 0) . . . ip + i(^ n \ 0)^+1 (^ n+1 ^, r) . . . i/j+iiC^, 0) over many trajectories (seedings) 
will effectively converge towards a discretized version of Eq. ( 1241) . Due to the formal equiva- 
lence of the semiclassical solutions to the quantum ones, the product ^_ 1 (^ 1 , r) n ^_ 1 (^ 2 , t) u 
will consequently converge to the quantum expectation value (|22|) . To find correlation func- 
tions of other operators, other seedings have to be used, which can be found in an analogous 
way. Table [I] summarizes these relationships. 

As far as densities and populations are concerned, the convergence of the averaging 
procedure towards the quantum solution is fairly fast. For instance, as depicted in Fig. |5l 
one typically has to average over a couple of thousand trajectories, to obtain well converged, 
smooth density profiles. We note that the curves lie generally below the quantum mechanical 
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Quantum Operator 


Seeded Fields 


MM 


^+i(£',o) 


^+i(M 


v^i(e',o) 
e+(oy) 


i-l(M 


^+i(e',o) 



TABLE I: Summary of relations between quantum mechanical expectation values and random 
initial seeds to the mean-field model. 

results. This is, however, a purely numerical effect, and it disappears when the number of 
trajectories increases. For a sufficiently large sample of trajectories, we hence also have true 
numerical convergence. Clearly, the averaging procedure is also applicable to correlations of 
higher order. One can calculate any normal ordered correlation function of the system by 
averaging over a sufficiently large ensemble of MF trajectories, provided that the operators 
involved in the correlation function have the same seeding requirements. The convergence, 
however, becomes rapidly slower with every order added, such that higher order correlations 
will require a larger number of trajectories. 



V. RESULTS BEYOND THE MEAN FIELD MODEL 



A. Population Ratio 

An easily accessible quantity in a BEC superradiance experiment is the ratio of backwards 
to forwards recoiling atoms. In the quantum model, the calculation of this quantity is 
straight forward and unambiguous. From Eqs. ( fT5l) and fl!2aft we find 

r/o Ad ^o d ^V(OI 2 ^(e)l 2 l^i(7^,r)pJ 

One cannot expect the MF model to deliver reliable results for such a quantity, due to its 
ambiguous initialization related to the seeding. To find out how different the predictions 
of the two models for this ratio are, the MF equations of motion were seeded according to 
Eq. (l20j) . and the predictions of both models for the time evolution of the ratio in the two 
regimes are plotted in Fig. [6j 
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FIG. 6: (Color online)Time evolution of the ratio of backwards to forwards recoiling atoms as a 
function of scaled time in the quantum (black, solid) and the MF (red, dashed) model, for the 
weak-pulse (a) and the strong pulse (b) regimes. 

Let us recall here that according to experimental observations, the population of the 
backwards atomic sidemodes is highly suppressed in the weak-pulse regime. Our quantum 
theoretical predictions reproduce these observations, i.e. the ratio is closer to one in the 
strong-pulse regime. In the MF model, it never exceeds 0.3 throughout its evolution in 
either of the two regimes, even though it does attain higher values for strong couplings. 
As discussed in Sec. [Til the suppression of backwards recoiling atoms due to the energy 
mismatch is expected to set in at r « r c . Indeed, as depicted in Fig. |6j in the weak 
pulse regime the ratio of backwards to forwards recoiling atoms drops for times r > r c . 
In the strong-pulse regime, however, and for the time scales consistent with the validity of 
our model, we are always well below r c , and the ratio increases monotonically. Nevertheless, 
even in this case, the onset of the suppression manifests itself in the temporal behavior of the 
growth rate of the ratio. Finally, it is also worth noting that the evolution of the ratio in the 
weak-pulse regime agrees qualitatively with the corresponding results in Q. A quantitative 
comparison, however, is rather difficult due to different definitions of the coupling strength. 
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B. Coherence of the scattered light 



Another class of questions typically addressed in a superradiance experiment pertains to 
the properties of emitted light. Of particular interest are the coherence properties which 
are described through correlation functions and are also accessible to measurements. For 
instance, the first-order correlation function G^\r,T) = (e + (A, r + T)e+(A, r)), describes 
the temporal coherence properties of the light that are relevant to an experiment of the 
Young's type, where the light at two times (i.e., at r and r + T) is superimposed to produce 
interference patterns . The visibility of the fringes in the observed pattern is proportional 
to the degree of first-order coherence, defined as 

V3(A,T + TP(A,r) 



On the other hand, intensity corre 
correlation function defined as 



ations are described through the normalized second order 



(% i rp\ (4 (A, r)e{(A, r + T)e + (A, r)e + (A, r + T))) 

This quantity is basically related to the probability of detecting a photon at time r+T, given 
that a photon was detected at time r. Definitions ( 12 7p and (j28p are general and applicable 
to any light source, but can be simplified considerably for stationary sources, where the 
properties of the light depend only on the delay time T. Actually, this is the case typically 
discussed in standard text books In the present setting, however, the process is by 

no means stationary, and thus well known expressions and conclusions are not necessarily 
applicable to our case. 

Using Eqs. ( lT5j) . one obtains 

^)( r , r + T ) = l + |^(r,r + T)| 2 , (29a) 

where 

G«(r,r + T)= / d^MO\ 2 Fi, (lA^r)F* ( lA ^r + T), (29b) 
Jo 

while t) is given by Eq. (|2~T|) . 

Equation (I29aj) is a typical property of so-called chaotic light sources jl^j], albeit in our 
case the source is non stationary. Indeed, as depicted in Fig. [71 even within the undepleted 
pump approximation adopted throughout this work, for a given delay time T, g^(r, r + T) 
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depends crucially on r. In view of Eq. fl2D, we have g^^r) = 2 for all r, which is a 
manifestation of intensity correlations. In other words, the endfire photons tend to appear 
bunched and thus detecting a photon at time r, significantly increases the probability of 
detecting another photon simultaneously On the other hand, as T — > oo we obtain g^> (r, r+ 
T) — > 1, indicating that the intensities are uncorrelated for large delay times. Typically, 
these asymptotic behaviors of g^ 2 \r,r + T) are defined with respect to the characteristic 
coherence time T c of the light under investigation (i.e., T — > and T — > oo refer to T < T c 
and T 3> T c , respectively). Unfortunately, the validity of the present model restricts us to 
relatively small times, and we cannot provide quantitative estimates of the coherence time 
T c . Nevertheless, we can still draw some conclusions about the behavior of g^ 2 \r, r + T) in 
the weak and in the strong pulse regimes. Before this, it is also worth pointing out that 
measuring the correlation function in dependence of the delay time T (irrespective of r) 
would facilitate any experiment considerably. In practice, this can be achieved, for instance, 
by forming blocks of data pertaining to various r but the same delay time T, and then 
estimate g^ 2 \T) based on these blocks. Formally, this procedure corresponds to the time 
averaged degree of second order coherence given by 

~® m = Jo" dr(4(A, r)4(A, r + T)e + (A, r)e + (A, r + T)) 

9 1 ] r o To drJ(A,r)a(A,r + T) ' 1 ] 

where r denotes the time over which experimental data are collected. This expression is 



analogous to volume integrated correlation functions used in 15|. Note that when there is no 
dependence on r, Eq. (!30j) reduces to the the standard expression of g^ 2 \T) for stationary 
sources fl^ |. 

As depicted in Fig. [8], the behavior of g^ 2 \T) is substantially different in the weak and 
the strong pulse regimes. While in the weak pulse regime the correlation function seems to 
decay slowly but steadily, in the strong pulse regime we clearly have two stages. The initial 
transient regime is characterized by a rapid drop of g( 2 \T), which is followed by a regime of 
very slow decay. To a good approximation g^ 2 \T) decreases linearly with increasing delay 
times, in both regimes. Moreover, according to Fig. El the tendency of photons to arrive 
in bunches is much lower in the strong pulse regime than in the weak pulse regime. This 
behavior can be attributed to the production of backwards recoiling atoms at the expense of 
endfire photons (see red curve). Intensity correlation function can thus viewed as a measure 
of the contribution of backwards recoiling atoms to the total number of scattering events. 
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FIG. 7: (Color online) Behavior of the photon intensity correlation (r, r + T) as a function of r 
for various fixed delays T as given in the figures. The dashed line in the top figure gives the same 
function neglecting the (—1) mode. 

In this context, we can also interpret the behavior of g( 2 \r,T + T) as a function of r. 
As depicted in Fig. [7J for both regimes there seems to be a systematic temporal behavior 
of g {2 Kr,r + T) for a given delay time. The correlation function decreases for short times, 
while for larger times it increases (at least in the weak pulse regime). Such a behavior 
reflects changes in the statistical properties of the source, which can be also associated 
with the production of backward recoiling atoms. Indeed, as depicted in Fig. [TJ neglecting 
the backward mode (— ) in our equations of motion, one finds a monotonic behavior of 



(2), 



r, r + T) (in good qualitative agreement with predictions for conventional superradiance 
16|). For short times, photons that have been scattered into the endfire mode are consumed 
during the production of backward recoiling atoms which, as discussed in Sec. [TTJ can take 
place for times shorter than r c . This is also confirmed by the fact that according to Fig. 
[7J in the weak pulse regime for a given T, g( 2 \r, r + T) exhibits a minimum for times very 
close to r c . In the strong-pulse regime, although our time scales are always well below r c , the 
onset of suppression manifests itself in the minimum of the intensity correlation function. 
To complete the picture, it is important to note here that superradiant Rayleigh scattering 
off BECs basically involves the mixing of optical and matter waves, which is a nonlinear 
process. Thus, any changes in the densities and/or populations of the fields that are mixed, 
are expected to affect considerably the statistics of the scattered light. 
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FIG. 8: (Color online) Integrated correlation function g( 2 \T) as a function of the delay time T 
for the weak (top) and strong (bottom) pulse regime. Dashed line in the top figure is the same 
function neglecting backwards recoiling atoms. 

VI. CONCLUSIONS 

We have discussed the early stage of superradiant Rayleigh scattering off atomic conden- 
sates, in the framework of a quantum model that takes into account spatial effects. Exact 
analytic solutions to this model can be expressed in terms of integrals involving Bessel func- 
tions, and are substantially different from the corresponding semiclassical solutions that 
have been obtained previously in the context of a related mean-field model. Nevertheless, 
the predictions of the two models about density profiles and growth rates are in reasonable 
qualitative agreement. An exception to this behavior is a strong suppression of photonic 
endfire modes at early times, which is underestimated in the mean-field treatment. For a 
large collection of condensed atoms, the effect of quantum fluctuations on various observ- 
able quantities can be obtained by averaging over ensembles of semiclassical trajectories. 
Each trajectory corresponds to the solution of the mean-field equations of motion, where 
an appropriate random seeding has been used. Hence histograms of a particular observable 
will reveal its distribution according to the quantum model. This technique will be used 
in future work to study photon delay time statistics. The quantum predictions for the ra- 
tio of backward to forward recoiling atoms is in qualitative agreement with experimental 
observations as well as other theoretical treatments. Finally, the present model enabled us 
to calculate the statistical behavior of scattered photons, which is qualitatively different in 
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the Kapitza-Dirac and Bragg regimes. This difference can be attributed to the suppression 
of photonic endfire modes due to the large number of back-scattered atoms at early times. 
In both regimes, the presence of backscattering distinguishes the photon statistics from the 
ones observed in conventional superradiance. 
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Appendix A: Solving the Equations of Motion 



To find solutions to Eqs. we first apply the Laplace transform with respect to r — > p 
to all three equations and find 

1 



Cr^pW+ifor)} =- [^i(e,0) + i^(O^ P {e+(e,r)} 

1 



p + 2i 



o) - i«Vo(0£r-*{s+& r )} 



(Ala) 
(Alb) 







X-Q7^r^p{e + (^, t)} = - pC T ^ p {e + (^ r)} + e+(f , 0) 



IK 



MS)^{$+i{M} + MZ)^p{$U{M} . (Aic 



After inserting (lAlal) and ( lAlbl) in Eq. (lAlcj) we are left with a differential equation of the 
form 

(A2) 



— £ r ^p{e + (^,r)} = -a(t,p)£ T ^p{e + (£,r)} - b(£,p), 
which has the solution 

We + (£,r)} = e-^«'rf^{e + (0 1 r)} - e'^^P) f ^ *) h{ ^ p ). 

Jo 

(A3) 

This expression can be inserted in Eqs. (lAlafAlbj) . at which point we have closed expressions 
for all three Laplace transformed fields. The remaining inversion of the Laplace transforms 
can be found using elementary techniques given in e.g. 17|. The solutions are expressed in 
terms of the following functions 
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F 0t0 (y, z) =5(z) + 9(^)y|/i (2y/yz) - Q{z)e~ 2iz ^-J x fey/yt 

zi I 1 (2^/y(z-z')"jJ 1 (2 y /yz) 



Q(z)y / dz'e 
Jo 



y/(z - z')z' 



F lt0 {y, z) =/ (2VP) 6(z) - e(z)^ / dz'^Ll \l^y{z-z> 

Jo yz' l 



Ji 2^' 



(A4a) 

(A4b) 
(A4c) 



F 0>1 (y, z) =e-**J (2y/yz) Q(z) + Q(z)y/y / dz J -==I 1 



j e 



-i2z' 



F ltl (y,z)=Q{z) / d^e-^'Jo 2^R 



z — Z' 



yjz- z 

Jo(2Vy^ 



z — z 



Joi^Vyz 1 



(A4d) 
(A4e) 



F 2fi (y,z)=J-h(2^)Q(z)-Q(z) / dz'e 



Z — Z 



F^(y,z)=e~ 2l \l-J 1 (2^)Q(z) + e(z) / dz'e 
V Jo 



/„-2iz' 



—l l \2^y{z-z')\.h^^% 

(A4f) 

yh^yiz-z')}^^}, 

(A4g) 



z — z' 



with Jj and Jj the zth Bessel function of the first kind and the ith modified Bessel function 
respectively. They satisfy the recursion relations 



d 

—F^ u (u,v) =F^ +1)U (u,v) - F^ u+l (u,v) 



d_ 

dv 



Ffj, tU (u,v) =2F^ hu (u,v) 



(A5) 



2iFaJu, v) =F u , v -i{u, v) - Fu-iJu, v) 



which are a consequence of properties of the Laplace transform. For the solutions with 
neglected backwards scattering, we use the functions 

(A6) 
(A7) 

F 2 (y,z) = 9(^/^(2 VP). (A8) 



F (y,z) = 5(z)+e(z)^h(2y/yz) 
F 1 (y,z) = e(z)J (2Vyi) 
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